# Replication Script for Note on Besley, Sturm and Persson (2010) #
################ Adi Dasgupta and Sean Ingham #####################

################## Set Working Directory ##########################
###################################################################

setwd <- choose.dir() ## Files will be loaded from and saved to the selected folder

################# Libraries and Functions #########################
###################################################################

library(foreign) ## To load .dta STATA dataset
library(sandwich) ## Clustered standard errors
library(lmtest) ## Clustered standard errors

VC <- function(fm, dfcw, cluster){ ## Function, by Mahmood Arai (see reference below), to estimate variance-covariance matrix adjusted for clustering of errors
	M <- length(unique(cluster)) ## see http://people.su.se/~ma/clustering.pdf, section 3
	N <- length(cluster)
	dfc <- (M/(M-1))*((N-1)/(N-fm$rank))
	u <- apply(estfun(fm),2,
	function(x) tapply(x, cluster, sum))
	vcovCL <- dfc*sandwich(fm, meat=crossprod(u)/N)*dfcw
	return(vcovCL) 
}

############### Load Main Dataset #################################
###################################################################

data <- read.dta("formatteddata.dta") ## Load main dataset 

############### Fig. 1 ############################################
###################################################################

### Three OLS models with taxes as the outcome variable

predictors.taxes.comp <- paste(c("compnorm","gip","demcontrol","repcontrol",
	"as.factor(state)",names(data)[181:231], names(data)[57:107]),
	sep="",collapse="+") ## Creating explanatory variable string for model with just political competition
form.taxes.comp <- as.formula(paste("share_taxes_inc", predictors.taxes.comp, 
	sep = " ~ ")) ## Creating model string
ols.taxes.comp <- lm(form.taxes.comp,data=data) ## OLS regression
ols.taxes.comp$vcov <- VC(ols.taxes.comp,1,model.frame(ols.taxes.comp)[,6]) ## Variance-covariance matrix adjusted for clustering

predictors.taxes.dem <- paste(c("normdem","gip","demcontrol","repcontrol",
	"as.factor(state)",names(data)[181:231], names(data)[57:107]),
	sep="",collapse="+") ## Creating explanatory variable string for model with just Democratic vote share
form.taxes.dem <- as.formula(paste("share_taxes_inc", predictors.taxes.dem, 
	sep = " ~ ")) ## Creating model string
ols.taxes.dem <- lm(form.taxes.dem,data=data) ## OLS regression
ols.taxes.dem$vcov <- VC(ols.taxes.dem,1,model.frame(ols.taxes.dem)[,6]) ## Variance-covariance matrix adjusted for clustering

predictors.taxes.both <- paste(c("normdem","compnorm","gip","demcontrol","repcontrol",
	"as.factor(state)",names(data)[181:231], names(data)[57:107]),
	sep="",collapse="+") ## Creating explanatory variable string for model with just Democratic vote share
form.taxes.both <- as.formula(paste("share_taxes_inc", predictors.taxes.both, 
	sep = " ~ ")) ## Creating model string
ols.taxes.both <- lm(form.taxes.both,data=data) ## OLS regression
ols.taxes.both$vcov <- VC(ols.taxes.both,1,model.frame(ols.taxes.both)[,7]) ## Variance-covariance matrix adjusted for clustering

#### Estimated marginal effect of Democratic vote share in these three models

maintax <- c(coefficients(ols.taxes.comp)[2],-coefficients(ols.taxes.comp)[2]) ## With only political competition included
normtax <- c(coefficients(ols.taxes.dem)[2],coefficients(ols.taxes.dem)[2]) ## With only Democratic vote share included
besttax <- c(coefficients(ols.taxes.both)[2] + coefficients(ols.taxes.both)[3],
	coefficients(ols.taxes.both)[2] - coefficients(ols.taxes.both)[3]) ## With both included

#### Standard errors of the estimated marginal effect of Democratic vote share

maintaxse <- c(sqrt(ols.taxes.comp$vcov[2,2]), sqrt(ols.taxes.comp$vcov[2,2]))
normtaxse <- c(sqrt(ols.taxes.dem$vcov[2,2]), sqrt(ols.taxes.dem$vcov[2,2]))
besttaxse <- c(sqrt(ols.taxes.both$vcov[2,2] + ols.taxes.both$vcov[3,3] + 2*ols.taxes.both$vcov[2,3]),
	sqrt(ols.taxes.both$vcov[2,2] + ols.taxes.both$vcov[3,3] - 2*ols.taxes.both$vcov[2,3]))

#### Marginal effect plots (figures 1, 2, and 3 in Dasgupta and Ingham)

coef.vec <- c(maintax[1],maintax[2],normtax[1],normtax[2],besttax[1],besttax[2]) ## Creating a vector of marginal effect estimates
se.vec <- c(maintaxse[1],maintaxse[2],normtaxse[1],normtaxse[2],besttaxse[1],besttaxse[2]) ## Creating a vector of standard error estimates
index <- seq(1,6,by=2) ## Estimates for Democratic vote share < .5

coef.vecR <- coef.vec[index] 
se.vecR <- se.vec[index]
coef.vecD <- coef.vec[-index]
se.vecD <- se.vec[-index]

x.axis <- c(1:length(coef.vecR))
x.axis2 <- c((1-.5):(length(coef.vecR)-.5))

m <- floor(-4*sd(data$share_taxes_inc,na.rm=T)); M <- ceiling(4*sd(data$share_taxes_inc,na.rm=T))

pdf(file="taxes.pdf",width=10,height=7)
plot(x.axis2, coef.vecR, type = "p", axes = F, xlab = "", ylab = "", pch = 19, cex = 1.2,
    xaxs = "r",ylim=c(m,M),col='red',xlim=c(.4,3.2),
	main="Marginal Effect of Democratic Party Vote Share on Taxes")
points(x.axis,coef.vecD,type='p',col='blue',pch=17, cex = 1.2)
abline(v=c(2.25,1.25))
abline(h=0,lty=2)
segments(x.axis2, coef.vecR-qnorm(.975)*se.vecR, 
	x.axis2, coef.vecR+qnorm(.975)*se.vecR, lwd =  1.1,lty=1)
segments(x.axis2 -.1,coef.vecR-qnorm(.95)*se.vecR, 
	x.axis2 +.1, coef.vecR-qnorm(.95)*se.vecR, lwd = 1.1,lty=1)
segments(x.axis2 -.1, coef.vecR+qnorm(.95)*se.vecR,
	x.axis2 +.1, coef.vecR+qnorm(.95)*se.vecR, lwd = 1.1,lty=1)
segments(x.axis, coef.vecD-qnorm(.975)*se.vecD,
	x.axis, coef.vecD+qnorm(.975)*se.vecD, lwd =  1.1,lty=1)
segments(x.axis -.1, coef.vecD-qnorm(.95)*se.vecD,
	x.axis +.1, coef.vecD-qnorm(.95)*se.vecD, lwd = 1.1,lty=1)
segments(x.axis -.1, coef.vecD+qnorm(.95)*se.vecD,
	x.axis +.1, coef.vecD+qnorm(.95)*se.vecD, lwd = 1.1,lty=1)
axis(2, at = seq(m,M,by=1), labels =  seq(m,M,by=1), tick = T, 
    cex.axis = .8)
box(bty = "o")
axis(1, at = c(2.75,1.75,.75), label = c(expression(paste(d[s][t], " : Yes, ", kappa[s][t], " : Yes")),
	expression(paste(d[s][t], " : Yes, ", kappa[s][t], " : No")),
	expression(paste(d[s][t], " : No, ", kappa[s][t], " : Yes"))), 
	tick = T, ,mgp = c(2,.6,0),
	cex.axis = 1.2) 
mtext("Model Specification", side=1, line=2.5)
mtext("Marginal Effect on Taxes as a Percentage of State Income", side=2, line=3)
legend(x="top", legend=c(expression(paste("Effect of ", d[s][t], " from 0 to.5")), 
	expression(paste("Effect of ", d[s][t], " from .5 to 1"))),
       col=c("red", "blue"), pch=c(19,17),lty=c(1,1), bty="n") 
text(x=c(x.axis2, x.axis), 
	y=c(coef.vecR+qnorm(.975)*se.vecR+.5, coef.vecD+qnorm(.975)*se.vecD+.5), 
	c(rep("[0,.5)",times=3),
	rep("(.5,1]",times=3)))
dev.off()

############### Fig. 2 ############################################
###################################################################

### Three OLS models with infrastructures spending as the outcome variable

predictors.inf.comp <- paste(c("compnorm","gip","demcontrol","repcontrol",
	"as.factor(state)",names(data)[181:231], names(data)[57:107]),
	sep="",collapse="+") ## Creating explanatory variable string for model with just political competition
form.inf.comp <- as.formula(paste("share_cap_exp", predictors.inf.comp, 
	sep = " ~ ")) ## Creating model string
ols.inf.comp <- lm(form.inf.comp,data=data) ## OLS regression
ols.inf.comp$vcov <- VC(ols.inf.comp,1,model.frame(ols.inf.comp)[,6]) ## Variance-covariance matrix adjusted for clustering

predictors.inf.dem <- paste(c("normdem","gip","demcontrol","repcontrol",
	"as.factor(state)",names(data)[181:231], names(data)[57:107]),
	sep="",collapse="+") ## Creating explanatory variable string for model with just Democratic vote share
form.inf.dem <- as.formula(paste("share_cap_exp", predictors.inf.dem, 
	sep = " ~ ")) ## Creating model string
ols.inf.dem <- lm(form.inf.dem,data=data) ## OLS regression
ols.inf.dem$vcov <- VC(ols.inf.dem,1,model.frame(ols.inf.dem)[,6]) ## Variance-covariance matrix adjusted for clustering

predictors.inf.both <- paste(c("normdem","compnorm","gip","demcontrol","repcontrol",
	"as.factor(state)",names(data)[181:231], names(data)[57:107]),
	sep="",collapse="+") ## Creating explanatory variable string for model with just Democratic vote share
form.inf.both <- as.formula(paste("share_cap_exp", predictors.inf.both, 
	sep = " ~ ")) ## Creating model string
ols.inf.both <- lm(form.inf.both,data=data) ## OLS regression
ols.inf.both$vcov <- VC(ols.inf.both,1,model.frame(ols.inf.both)[,7]) ## Variance-covariance matrix adjusted for clustering

#### Estimated marginal effect of Democratic vote share in the three models

maininf <- c(coefficients(ols.inf.comp)[2],-coefficients(ols.inf.comp)[2]) ## With only political competition included
norminf <- c(coefficients(ols.inf.dem)[2],coefficients(ols.inf.dem)[2]) ## With only Democratic vote share included
bestinf <- c(coefficients(ols.inf.both)[2] + coefficients(ols.inf.both)[3],
	coefficients(ols.taxes.both)[2] - coefficients(ols.inf.both)[3]) ## With both included

#### Standard errors of the estimated marginal effect of Democratic vote share

maininfse <- c(sqrt(ols.inf.comp$vcov[2,2]), sqrt(ols.inf.comp$vcov[2,2]))
norminfse <- c(sqrt(ols.inf.dem$vcov[2,2]), sqrt(ols.inf.dem$vcov[2,2]))
bestinfse <- c(sqrt(ols.inf.both$vcov[2,2] + ols.inf.both$vcov[3,3] + 2*ols.inf.both$vcov[2,3]),
	sqrt(ols.inf.both$vcov[2,2] + ols.inf.both$vcov[3,3] - 2*ols.inf.both$vcov[2,3]))

#### Marginal effect plot

coef.vec <- c(maininf[1],maininf[2],norminf[1],norminf[2],bestinf[1],bestinf[2]) ## Creating a vector of marginal effect estimates
se.vec <- c(maininfse[1],maininfse[2],norminfse[1],norminfse[2],bestinfse[1],bestinfse[2]) ## Creating a vector of standard error estimates
index <- seq(1,6,by=2) ## Estimates for Democratic vote share < .5

coef.vecR <- coef.vec[index] 
se.vecR <- se.vec[index]
coef.vecD <- coef.vec[-index]
se.vecD <- se.vec[-index]

x.axis <- c(1:length(coef.vecR))
x.axis2 <- c((1-.5):(length(coef.vecR)-.5))

m <- floor(-2*sd(data$share_cap_exp,na.rm=T)); M <- ceiling(2*sd(data$share_cap_exp,na.rm=T))

pdf(file="infrastructure.pdf",width=10,height=7)
plot(x.axis2, coef.vecR, type = "p", axes = F, xlab = "", ylab = "", pch = 19, cex = 1.2,
    xaxs = "r",ylim=c(m,M),col='red',xlim=c(.4,3.2),
	main="Marginal Effect of Democratic Party Vote Share on Infrastructure Spending")
points(x.axis,coef.vecD,type='p',col='blue',pch=17, cex = 1.2)
abline(v=c(2.25,1.25))
abline(h=0,lty=2)
segments(x.axis2, coef.vecR-qnorm(.975)*se.vecR, 
	x.axis2, coef.vecR+qnorm(.975)*se.vecR, lwd =  1.1,lty=1)
segments(x.axis2 -.1,coef.vecR-qnorm(.95)*se.vecR, 
	x.axis2 +.1, coef.vecR-qnorm(.95)*se.vecR, lwd = 1.1,lty=1)
segments(x.axis2 -.1, coef.vecR+qnorm(.95)*se.vecR,
	x.axis2 +.1, coef.vecR+qnorm(.95)*se.vecR, lwd = 1.1,lty=1)
segments(x.axis, coef.vecD-qnorm(.975)*se.vecD,
	x.axis, coef.vecD+qnorm(.975)*se.vecD, lwd =  1.1,lty=1)
segments(x.axis -.1, coef.vecD-qnorm(.95)*se.vecD,
	x.axis +.1, coef.vecD-qnorm(.95)*se.vecD, lwd = 1.1,lty=1)
segments(x.axis -.1, coef.vecD+qnorm(.95)*se.vecD,
	x.axis +.1, coef.vecD+qnorm(.95)*se.vecD, lwd = 1.1,lty=1)
axis(2, at = seq(m,M,by=1), labels =  seq(m,M,by=1), tick = T, 
    cex.axis = .8)
box(bty = "o")
axis(1, at = c(2.75,1.75,.75), label = c(expression(paste(d[s][t], " : Yes, ", kappa[s][t], " : Yes")),
	expression(paste(d[s][t], " : Yes, ", kappa[s][t], " : No")),
	expression(paste(d[s][t], " : No, ", kappa[s][t], " : Yes"))), 
	tick = T, ,mgp = c(2,.6,0),
	cex.axis = 1.2) 
mtext("Model Specification", side=1, line=2.5)
mtext("Marginal Effect on Infrastructure Spending as a Pecentage of the Budget", side=2, line=3)
legend(x="top", legend=c(expression(paste("Effect of ", d[s][t], " from 0 to.5")), 
	expression(paste("Effect of ", d[s][t], " from .5 to 1"))),
       col=c("red", "blue"), pch=c(19,17),lty=c(1,1), bty="n") 
text(x=c(x.axis2, x.axis), 
	y=c(coef.vecR+qnorm(.975)*se.vecR+1, coef.vecD+qnorm(.975)*se.vecD+1), 
	c(rep("[0,.5)",times=3),
	rep("(.5,1]",times=3)))
dev.off()

############### Fig. 3 ############################################
###################################################################

### Three OLS models with right-to-work laws as the outcome variable

predictors.rtw.comp <- paste(c("compnorm","gip","demcontrol","repcontrol",
	"as.factor(state)",names(data)[160:231], names(data)[36:107]),
	sep="",collapse="+") ## Creating explanatory variable string for model with just political competition
form.rtw.comp <- as.formula(paste("rtw", predictors.rtw.comp, 
	sep = " ~ ")) ## Creating model string
ols.rtw.comp <- lm(form.rtw.comp,data=data) ## OLS regression
ols.rtw.comp$vcov <- VC(ols.rtw.comp,1,model.frame(ols.rtw.comp)[,6]) ## Variance-covariance matrix adjusted for clustering

predictors.rtw.dem <- paste(c("normdem","gip","demcontrol","repcontrol",
	"as.factor(state)",names(data)[160:231], names(data)[36:107]),
	sep="",collapse="+") ## Creating explanatory variable string for model with just Democratic vote share
form.rtw.dem <- as.formula(paste("rtw", predictors.rtw.dem, 
	sep = " ~ ")) ## Creating model string
ols.rtw.dem <- lm(form.rtw.dem,data=data) ## OLS regression
ols.rtw.dem$vcov <- VC(ols.rtw.dem,1,model.frame(ols.rtw.dem)[,6]) ## Variance-covariance matrix adjusted for clustering

predictors.rtw.both <- paste(c("normdem","compnorm","gip","demcontrol","repcontrol",
	"as.factor(state)",names(data)[160:231], names(data)[36:107]),
	sep="",collapse="+") ## Creating explanatory variable string for model with just Democratic vote share
form.rtw.both <- as.formula(paste("rtw", predictors.rtw.both, 
	sep = " ~ ")) ## Creating model string
ols.rtw.both <- lm(form.rtw.both,data=data) ## OLS regression
ols.rtw.both$vcov <- VC(ols.rtw.both,1,model.frame(ols.rtw.both)[,7]) ## Variance-covariance matrix adjusted for clustering

#### Estimated marginal effect of Democratic vote share in the three models

mainrtw <- c(coefficients(ols.rtw.comp)[2],-coefficients(ols.rtw.comp)[2]) ## With only political competition included
normrtw <- c(coefficients(ols.rtw.dem)[2],coefficients(ols.rtw.dem)[2]) ## With only Democratic vote share included
bestrtw <- c(coefficients(ols.rtw.both)[2] + coefficients(ols.rtw.both)[3],
	coefficients(ols.rtw.both)[2] - coefficients(ols.rtw.both)[3]) ## With both included

#### Standard errors of the estimated marginal effect of Democratic vote share

mainrtwse <- c(sqrt(ols.rtw.comp$vcov[2,2]), sqrt(ols.rtw.comp$vcov[2,2]))
normrtwse <- c(sqrt(ols.rtw.dem$vcov[2,2]), sqrt(ols.rtw.dem$vcov[2,2]))
bestrtwse <- c(sqrt(ols.rtw.both$vcov[2,2] + ols.rtw.both$vcov[3,3] + 2*ols.rtw.both$vcov[2,3]),
	sqrt(ols.rtw.both$vcov[2,2] + ols.rtw.both$vcov[3,3] - 2*ols.rtw.both$vcov[2,3]))

#### Marginal effect plot

coef.vec <- c(mainrtw[1],mainrtw[2],normrtw[1],normrtw[2],bestrtw[1],bestrtw[2]) ## Creating a vector of marginal effect estimates
se.vec <- c(mainrtwse[1],mainrtwse[2],normrtwse[1],normrtwse[2],bestrtwse[1],bestrtwse[2]) ## Creating a vector of standard error estimates
index <- seq(1,6,by=2) ## Estimates for Democratic vote share < .5

coef.vecR <- coef.vec[index] 
se.vecR <- se.vec[index]
coef.vecD <- coef.vec[-index]
se.vecD <- se.vec[-index]

x.axis <- c(1:length(coef.vecR))
x.axis2 <- c((1-.5):(length(coef.vecR)-.5))

m <- floor(-4*sd(data$rtw,na.rm=T)); M <- ceiling(4*sd(data$rtw,na.rm=T))

pdf(file="rtw.pdf",width=10,height=7)
plot(x.axis2, coef.vecR, type = "p", axes = F, xlab = "", ylab = "", pch = 19, cex = 1.2,
    xaxs = "r",ylim=c(m,M),col='red',xlim=c(.4,3.2),
	main="Marginal Effect of Democratic Party Vote Share on Right-to-work Laws")
points(x.axis,coef.vecD,type='p',col='blue',pch=17, cex = 1.2)
abline(v=c(2.25,1.25))
abline(h=0,lty=2)
segments(x.axis2, coef.vecR-qnorm(.975)*se.vecR, 
	x.axis2, coef.vecR+qnorm(.975)*se.vecR, lwd =  1.1,lty=1)
segments(x.axis2 -.1,coef.vecR-qnorm(.95)*se.vecR, 
	x.axis2 +.1, coef.vecR-qnorm(.95)*se.vecR, lwd = 1.1,lty=1)
segments(x.axis2 -.1, coef.vecR+qnorm(.95)*se.vecR,
	x.axis2 +.1, coef.vecR+qnorm(.95)*se.vecR, lwd = 1.1,lty=1)
segments(x.axis, coef.vecD-qnorm(.975)*se.vecD,
	x.axis, coef.vecD+qnorm(.975)*se.vecD, lwd =  1.1,lty=1)
segments(x.axis -.1, coef.vecD-qnorm(.95)*se.vecD,
	x.axis +.1, coef.vecD-qnorm(.95)*se.vecD, lwd = 1.1,lty=1)
segments(x.axis -.1, coef.vecD+qnorm(.95)*se.vecD,
	x.axis +.1, coef.vecD+qnorm(.95)*se.vecD, lwd = 1.1,lty=1)
axis(2, at = seq(m,M,by=1), labels =  seq(m,M,by=1), tick = T, 
    cex.axis = .8)
box(bty = "o")
axis(1, at = c(2.75,1.75,.75), label = c(expression(paste(d[s][t], " : Yes, ", kappa[s][t], " : Yes")),
	expression(paste(d[s][t], " : Yes, ", kappa[s][t], " : No")),
	expression(paste(d[s][t], " : No, ", kappa[s][t], " : Yes"))), 
	tick = T, ,mgp = c(2,.6,0),
	cex.axis = 1.2) 
mtext("Model Specification", side=1, line=2.5)
mtext("Marginal Effect on Existence of Right-to-work Laws", side=2, line=3)
legend(x="top", legend=c(expression(paste("Effect of ", d[s][t], " from 0 to.5")), 
	expression(paste("Effect of ", d[s][t], " from .5 to 1"))),
       col=c("red", "blue"), pch=c(19,17),lty=c(1,1), bty="n") 
text(x=c(x.axis2, x.axis), 
	y=c(coef.vecR+qnorm(.975)*se.vecR+.2, coef.vecD+qnorm(.975)*se.vecD+.2), 
	c(rep("[0,.5)",times=3),
	rep("(.5,1]",times=3)))
dev.off()